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Abstract We study the dynamic behaviour of ID spring- 
block models of friction when the external loading is ap- 
plied from a side, and not on all blocks like in the classical 
Burridge-Knopoff-like models. Such a change in the loading 
yields specific difficulties, both from numerical and phys- 
ical viewpoints. To address some of these difficulties and 
clarify the precise role of a series of model parameters, we 
start with the minimalistic model by Maegawa et al. (Tri- 
bol. Lett. 38 313, 2010) which was proposed to reproduce 
their experiments about precursors to frictional sliding in 
the stick-slip regime. By successively adding (i) an internal 
viscosity, (ii) an interfacial stiffness and (iii) an initial tan- 
gential force distribution at the interface, we manage to (i) 
avoid the model's unphysical stress fluctuations, (ii) avoid 
its unphysical dependence on the spatial resolution and (iii) 
improve its agreement with the experimental results, respec- 
tively. Based on the behaviour of this improved ID model, 
we develop an analytical prediction for the length of pre- 
cursors as a function of the applied tangential load. We also 
discuss the relationship between the microscopic and macro- 
scopic friction coefficients in the model. 
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1 Introduction 

The dynamics of frictional interfaces are crucial to many sit- 
uations in mechanical engineering geosciences |3J] or 
biology |0, Etl- Today, after decades of studies, the science of 
contacts under time-invariant loading conditions, e.g. static 
contacts or steady sliding contacts [6-8], has reached a high 
level of advancement which, in manyinstances, enables quan- 
titative reproduction of global M or local 0, [13, LlH] 
measurements. In contrast, the dynamics of contacts under 
rapidly evolving loads or during fast unstable motion like 
stick-slip Ill2l4l4ll is far less understood. In particular, recent 
experiments on the transition from static to kinetic friction 
of side-driven poly(methyl methacrylate) (PMMA) rough 
samples forming line contacts with a rough PMMA sub- 
strate have revealed unexpected features Ill5l4l8n . The tran- 



sition occurs through the fast (comparable to the speed of 
soun d) prop agation of micro-slip fronts through the con- 
tact 11151 llal . It can also be preceded by a series of fronts 
that arrest before havi ng rup tured the whole contact, thus 
denoted as precursors | T3 12 1 • 

These results, which may have important implications 
for e.g. the study of earthquakes, have triggered an active 
modelling activity. Braun et al. Jl9ll . using a one-dimensional 
(ID) spring -block model with a complex time-dependent 
friction law, produced three types of micro-slip front veloci- 
ties, analogous to that observed in 1 15 ] . Maegawa et al. lfl7ll , 
using a ID spring -block model with a simple Amontons- 
Coulomb (AC) friction law, showed that the length of pre- 
cursors is modified when the external normal load is made 
asymmetric. Scheibert and Dysthe J20I1 . using a quasi-static 
ID model with AC friction, showed how the increasing tan- 
gential load itself induces an increasing pressure asymmetry 
which influences the precursors' series. Due to the intrin- 
sic limitations of ID models, all these studies yielded only a 
qualitative agreement with experiments. Very recently, Tr0mbor£ 
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et al. demonstrated, using a 2D spring-block model with 
AC friction, that quantitative agreement with the kinemat- 
ics {i.e. the properties of the states in which no micro-slip 
front is propagating) of the experiments requires an accu- 
rate description of the interfacial stresses, and therefore the 
use of realistic boundary conditions on the sample 1121 1. 

Compared to such a 2D model, the strongest advantage 
of ID models is that their results are much easier to analyse 
and understand, so that they provide opportunities for the- 
oretical approaches. They also require simpler implementa- 
tion and lower computational power. In most situations, in 
which very accurate results are not needed and/or a thor- 
ough qualitative understanding of the behaviour of the sys- 
tem is desired, ID models are preferable. From the pioneer- 
ing work of Burridge and Knopoff |22j], spring -block mod- 



els of friction have been extensively studied (see e.g. 1123 
l26ll ). These models have mainly been used to describe the 
statistical properties of the series of earthquakes at a seis- 
mic fault. Inertial blocks are connected in series via internal 
springs that model the crust's elasticity. A homogeneous tec- 
tonic loading is modelled by coupling, via loading springs, 
all individual blocks to the same rigid driving body. Such 
statistical analysis of homogeneously driven systems con- 
trast with the recent ID studies by Braun et al. i\% and Mae- 
gawa et al. lfl7ll , in which the time evolution of side-driven 
systems is analysed deterministic ally, in order to produce 
data comparable to the experimental measurements. The two 
models by Braun et al. and Maegawa et al. are actually 
very different. Both consider an array of blocks connected 
by internal springs, but the model by Braun et al. Il9ll also 
considers viscous dissipation and a complex time-dependent 
friction law emerging from the collective behaviour of in- 
terfacial springs with random stiffness, breaking threshold 
and reattachment time. In contrast the model by Maegawa 
et al. IT7I] . which only considers blocks and springs and the 
minimalistic AC friction law, is probably the simplest possi- 
ble model. 

On the one hand, we will see that, due to its extreme 
simplicity, the model by Maegawa et al. fnll yields results 
that are strongly resolution dependent, which prevents ro- 
bust comparison with experiments. On the other hand, it is 
difficult to disentangle the respective roles of the many pa- 
rameters of the model by Braun et al. 11911 . The scope of this 
article is therefore to (i) construct, step by step, a minimal 
ID side-driven spring-block model, the results of which are 
essentially resolution independent and (ii) qualify the capa- 
bilities of this model to reproduce the main qualitative fea- 
tures of recent experimental observations. 

This Letter is organised as follows. In section 2, we de- 
scribe the model by Maegawa et al. lfl7ll and show its limi- 
tations. In section 3, we improve this model by introducing 
successively an internal viscosity, an interfacial stiffness and 
an initial tangential force distribution. In section 4, we de- 



scribe an analytic prediction for the length of precursors in 
both Maegawa et al.'s and our improved model. 



2 The model by Maegawa et al. and its limitations 



2. 1 The model by Maegawa et al. 11711 



In the model developed by Maegawa et al. 11711 . the slider is 
modelled as a chain of blocks connected by springs (Fig.Q]), 
with material spring constant k and block mass m = M/N, 
where M is the total mass of the slider and N is the num- 
ber of blocks. In experiments, the base (also called track) is 
fixed on a very stiff support, and it is therefore modelled as a 
rigid surface for simplicity. The tangential force is applied at 
the trailing edge of the system through a loading spring with 
stiffness K. One end of this spring is attached to the trail- 
ing edge block (block 1), while the other end of the spring 
moves at a constant velocity V. The normal force pi is im- 
posed on each block, satisfying the criterion Ytf=i p\ = F^, 
where F^ is the total applied normal force. 
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Fig. 1 Schematics of the model system. The sample is modeled by N 
blocks of mass m connected in series through springs of stiffness k. The 
trailing edge of the system (block 1) is slowly driven through a loading 
spring of stiffness K. Each block is also submitted to a normal force p, 
and to a friction force /,. 



The equations of motion are given by 

{k{u2 — «i) +Fr + fi, 7i = l 
k(u n+ \ - 2u„ + m„_i) +/„, 2 < n < N— 1 (1) 
k(u N -i - u N )+f N , n =N, 

where u n = u n (t) is the position of block n as a function 
of time relative to its equilibrium position and "denotes the 
double derivative with respect to t . Fj = Fj{f) is the driving 
force (or tangential force/load) given by 



Ft =K(Vt -Mi). 



(2) 



A local friction law giving the friction forces /„ is im- 
posed between the blocks and the base. AC friction is used 
with local kinetic and static friction coefficients jU^ and jj, s , 
respectively. The resulting global friction coefficients are 
similarly denoted by )1k and /a^. The friction force on block 
«, /„, is then given by 



fn = 



< HsPn, 







-sgn(u n )n k p n , ii„ 7^0, 



(3) 
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where, when u„ = 0, equilibrium of block n imposes that /„ 
balances all other forces acting on block n. 

The material spring constant k is chosen such that the 
elastic deformation of the model is similar to that of a linear 
elastic medium with Young's modulus E, which yields: 



k=(N-l)ES/L, 



(4) 



where S and L are the cross-section area and the length of 
the slider, respectively. 



In the experiments by Maegawa et al. 11711 . an asymmet- 
ric normal loading was used, leading to the following linear 
model for p„: 



F N ( 2n-N-\ 

Pn = 1 9 

' N \ N-l 



(5) 



where E [—1, 1] is a measure of the non-uniformity in the 
normal loading. 

The values of the parameters are chosen to be in agree- 
ment with their respective values in the experiments by Mae- 
gawa et al. El, which are: K = 0.8 MN/m, V = 0.1 mm/s, 
F N = 400 N, M = 0.012 kg, L = 100 mm, S = 100 mm 2 , 
E = 2.5 GPa, n s = 0.7 and pt k = 0.45. 



2.2 Results and limitations of the model by Maegawa et al. 

We have implemented the model by Maegawa et al. and 
tested our code by comparing it to (i) two other similar codes 
and (ii) the analytical solution of the equation of motion for a 
one-block system. The simulation starts with each spring at 
its equilibrium length, i.e. u„ — for all blocks. From t = 0, 
the force from the driving spring on block 1 is increased. 
Figure [2^ shows that the resulting time evolution of Fr(t) 
exhibits stick-slip behaviour. Each time a drop in the driv- 
ing force is observed, some part of the slider moves relative 
to the base. We call the short time intervals during which 
movement occurs events. 

Between events, no block is moving. Since the driving 
force is applied only at the trailing edge block, only block 
1 is loaded and eventually reaches its static friction thresh- 
old, so that all events must nucleate at the trailing edge. The 
movement of block 1 then loads block 2, which itself reaches 
its threshold and so on. This succession of blocks starting to 
move defines a micro-slip front, which propagates towards 
the leading edg e, in analogy with the fronts observed in ex- 
periments lll5l Il7ll . The distance from the trailing edge to 



the micro-slip front as a function of time, x/(t ), is shown in 
Fig- Eh- If this front reaches the leading edge, the event is a 
global event, and the whole slider moves relative to the base. 

From Fig. [2] it is evident that not all events are global: 
smaller events are observed between global events. In addi- 
tion, a series of events with increasing maxima of is 
seen to precede the first global event. These events occur for 
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Fig. 2 Time evolution of Fj (a) and xt (b) in the model by Maegawa et 
al. using N = 10 (as in fl7ll ) and = 0. Before macroscopic stick-slip 
(reached when t rs 2.8 s), the loading curve is punctuated by partial 
relaxations associated with precursors to sliding, i.e. micro-slip fronts 
spanning a length smaller than L. 



Ft well below the macro scop ic static friction threshold, and 
are called precursors lfl6ill7ll . The maximum of x/(t) during 
a precursor event, i.e. the length of a precursor, is denoted by 

To perform a quantitative comparison between their ex- 
perimental and numerical results, Maegawa et al. focused 
on the relationship between the normalized length L p /L of 
the series of precursors (Fig.|2j)) and the normalized tangen- 
tial force Fj/Fn at which they are triggered (Fig.|2^). To do 
this, they discarded all simulated precursors having a length 
smaller than any of the previous events, with the justification 
that no such smaller event was observed in the experiments. 
Note that 2D models produce series of precursors of mono- 
tonically increasing length, so that practically no event has 
to be discarded 112 111 . 

Figure [3] shows the tangential force T normalised by the 
normal force p on each block at the time of initiation and ar- 
rest of a precursor of length L p = Q.6L. The tangential force 
is here defined as the total force on a block excepting the 
friction force. Each event is initiated when the tangential 
force on block 1 reaches the local static friction threshold. 
As block 1 moves, the tangential force on block 2 increases, 
eventually reaching the local static friction threshold, and 
starts to move. The precursor event arrests when the tangen- 
tial force built on a block by its left neighbour is not suffi- 
cient to make it reach its static friction threshold. The slow 
loading of block 1 continues, and will eventually trigger a 
new event nucleating at the trailing edge. 

The spatial resolution used in the above results, i.e. N = 
10 as used by Maegawa et al. [17:1], is rather low. This is es- 
pecially evident in Fig. [3] which only contains 10 data points 
for x/p. We expected to improve these results by simply in- 
creasing the resolution to N = 100. The corresponding load- 
ing and front position curves, Fj{t) and Xf(t ), are shown in 
Fig. a 
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Fig. 3 Tangential force distribution at initiation of a microslip-event 
(o). Block 1 has reached its slip threshold (f/p = fa). The event in- 
volves all blocks from 1 to 6 (i.e. a precursor of length L„/L = 0.6), 
leading to a modified tangential force distribution at arrest (□). Results 
obtained using N = 10 (as in ITHl ) and = 0. 



Or 




Fig. 4 Time evolution of Fj (a) and xt (b) in the model by Maegawa et 
al. using N = 100 and = 0. Comparison with Fig. [2] shows a drastic 
increase in the time-frequency of micro-slip events. 



One of the limitations of the model by Maegawa et al. 
is now clearly evident. By only changing the spatial reso- 
lution, the loading and front position curves are changed 
significantly. Some aspects are unchanged: precursors pre- 
cede the first global event and then stick-slip behaviour is 
observed. The final average level of Fj / f/v also appears to 
be conserved. However, the amplitude of the drops in the 
loading curve is reduced, while the number of events, both 
global and precursory, is seen to increase significantly. 

To illustrate this scaling with respect to the model reso- 
lution, we have plotted in Fig. [5] the evolution with N of (i) 
the total number of events and the number of global events 
(Fig. |5^), and (ii) the total number of precursors and the 
number of precursors longer than any previous ones (Fig.|5p). 
An approximately linear increasing trend is observed in all 
four curves. This behaviour is problematic as soon as one 
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Fig. 5 Number of different kinds of events as a function of N, us- 
ing 0=0. (a) Number of global events (solid line) for f 6 [5 s,20 s] 
and total number of events (dashed line), (b) Number of precursors 
longer than any previous one (solid line) and total number of precur- 
sors (dashed line). All curves are stongly increasing functions of N, 
indicating unphysical resolution-dependence of the model's results. 



wants to compare to experiments, in which the size of the 
drops in Ft/Fn and the number of events are well-defined 
experimental measurements. Maegawa et al. |l7j| used N = 
10, which produced a number of precursors similar to that 
observed in their experiments, but this agreement appears 
to be casual. A robust model should produce almost identi- 
cal numbers of events whatever the spatial discretization of 
the slider. In this respect, note that 2D models do produce a 
resolution-independent number of events 12111 . 

Another problematic behaviour of the model affects the 
tangential force spatial distribution, as shown in Fig. [6] The 
tangential force has been plotted at three different times: at 
t = 0.5 s, t = 3 s, and at the arrest of the precursor of length 
Lp/L = 0.7. Strong oscillations are observed, with a half 
spatial period of the order of the lattice spacing, whatever 
the number of blocks N used. Again this unphysical depen- 
dence on the model resolution impedes a robust comparison 
with experimental measurements of the tangential stress dis- 
tribution t(x), like e.g. those of II 1 811 - 

In their quantitative comparison between model and ex- 
perimental results, Maegawa et al. IU7I1 focused on the length 
of precursors L p as a function of the tangential load Ft, 
L p (Ft). This relationship has also been studied experimen- 
tally by Rubinstein et al. JTfj]. While similar behaviours are 
found, we will only compare our model results to the exper- 
imental results of Maegawa et al.. 

The experimental setup by Maegawa et al. allowed for 
non-uniformities in the normal loading, and they studied its 
consequences on the length of precursors. The non-uniformity 
in the normal loading is modelled as an asymmetric distri- 
bution of the normal loads p„ by using Eq. (f5J). The value of 
the parameter 6 = ±0.833 is chosen to be in agreement with 
its corresponding experimental value. 

The length of precursors L p is plotted as a function of 
the tangential force Fj at event arrest in Fig. [7] Three dif- 
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Fig. 6 Tangential force distribution % n normalised by the normal force 
p n on each block, at three different times: at t = 0.5 s, t = 3 s, and at the 
arrest of the precursor of length L„/L = 0.7 using N = 100 and 9 = 
0. Unwanted large oscillations with a half-period equal to the lattice 
spacing appear. 




Fig. 7 Length of precursors L p normalised by the system length L as a 
function of the tangential load Fj at event arrest normalised by the nor- 
mal load F N for N = 10 and 9 = 0.833 (OX 6 = (□) and 9 = -0.833 
(0), and N = 100 and 9 = 0.833 (*), 9 = (x) and 9 = -0.833 (+). 
Solid lines are the analytical predictions described in Sec. [4] 



ferent values of 9 are used, and the results for both N = 10 
and N — 100 are included. As discussed above, the number 
of precursors increases with N. However, the shape of the 
curves does not change significantly with N, thus enabling 
comparison with the shape of the experimental curves. 

The qualitative behaviour of L p as a function of Ft is 
in agreement with experiments. For 9 = 0.833, the normal 
force on the trailing edge is reduced, leading to longer pre- 
cursors for the same tangential load compared to 0=0. 
On the other hand, 9 = —0.833 leads to an increased nor- 
mal load on the trailing edge, and precursors are shorter 
for a given Fj. All curves converge to the same point at 
L p /L= 1, meaning that the global static friction threshold 
jj-s (the value of Ft /Fn when the first global event occurs) 



is independent of the normal loading distribution. The value 
of jig appears to be approximately 0.45, which incidentally 
is the value of the local kinetic friction coefficient jU^. The 
reason for this will become clear in Sec. [4] A quantitative 
comparison between the experimental results in B17T1 and 
the model results in Fig. [7] however, reveals large discrepan- 
cies: all three experimental curves are found way below their 
simulated counterpart, meaning that the simulation strongly 
overestimates L p for any given Ft; the rapid increase in pre- 
cursor length after L p /L w 0.5 that is observed experimen- 
tally has no equivalent in the model. 

Summing up, three main limitations of the model by 
Maegawa et al. have been identified: (i) the tangential force 
shows large oscillations, the wavelength of which is con- 
trolled by the lattice spacing, (ii) the number of all kinds of 
events is an increasing function of N and (iii) the quantitative 
agreement with the L p (Ft) -curve between model and ex- 
periments is poor. In the following section, we propose im- 
provements of the model that contribute to overcome these 
limitations. 



3 Improvements of the model by Maegawa et al. 

3.1 Introducing a relative viscous damping 

Resolution dependent oscillations are known to occur in Burridge- 
Knopoff-like models and more generally in dynamic rupture 
models involving AC friction at the interface between dis- 
similar elastic media J27I1 . Classical ways to reduce them 
significantly are either to regularize the AC friction law (see 
e.g. 112811 and references therein) or to introduce a viscous 
damping in the system ll9U2lLl29l43lll . Both are physically 
sound, but we choose to adopt the second approach. Phys- 
ically, such viscosity is a way to model the energy dissipa- 
tion that any material undergoes during deformation. After 
Knopoff and Ni 1310 , we consider the following form for the 
viscous force F% : 



77 (ii 2 -«i) , n = 1 

F,? = < 1? (w«+i - 2m„ + u n -i) , 2 < n < N- 1 
77 (%_i - %) , n=N, 



(6) 



which is a damping on the relative motion of neighbouring 
blocks. As in e.g. 11171 12111 . we assume that energy dissipa- 
tion due to the motion of a block relative to the substrate is 
satisfactorily included in /„, and therefore do not, in contrast 
to e.g. §19, 32], include any viscous damping at the interface 
in our system. This also serves to keep the model as simple 
as possible. The equations of motion are then given by 

(k(u 2 -m)+F T +F?+f h n=\ 
mu„ = < k(u n+ 1 - 2u„ + u„- 1 ) + F,? + f„ , 2 < n < N - 1 
[ k(u N -\ - u N ) +F$ + fy, n = N. 

(7) 
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Fig. 8 The tangential force T normalised by the normal force p on 
each block is plotted at three different times: at t = 0.5 s, t = 3 s, and 
at the arrest of the precursor of length L p /L = 0.7 from a simulation 
including the relative viscous damping and using N = 100, = and 
7] = V0.1 \fkm. Comparison with Fig. [6] shows that lattice-controlled 
fluctuations have disappeared. 



The tangential force T„ is still defined as the sum of all 
forces on a block excepting the friction force, i.e. now in- 
cluding F„ . It remains to choose a reasonable value of the 
damping coefficient 77 . 

In App.|A] the value of 77 for critical damping of waves 
of wavelength A = 2a, corresponding to the cut-off wave- 
length, is calculated. The result is rj ( = \fkm~. It is also shown 
that the values of 77 corresponding to critical damping of 
higher wavelength oscillations are always larger than \fkm~. 
We want to damp out oscillations of wavelength A = 2a, 
and using 77 = \fkrn is then a possibility. However, waves 
of other wavelengths close to A = 2a will also be highly 
damped, causing significant changes to the dynamics. Since 
this is an undesirable effect, a compromise has to be made. 
As suggested in Jjlll . the value 



QAV km ^ 0.32V km 



(8) 



is used in the following. Note that since k^N- 1 and m °= 
1 /N, 77 is TV-independent for N » 1 . 

Figure[8]shows the tangential force distribution obtained 
when relative viscous damping is included. The improve- 
ment with respect to Fig. [6] is clear: the short wavelength 
oscillations have almost disappeared, resulting in a physi- 
cally reasonable smooth tangential force profile, xjp also 
appears to be on average approximately equal to in rup- 
tured regions, a fact that will be utilised below to predict the 
precursor lengths. 

Some small one- or two-node spikes remain. They have 
a different origin, as they are mainly caused by the discrete- 
ness of the friction law: at the tip of the rupture, one block is 
moving and therefore increasing the force on its neighbour, 
which is still stationary. The region beyond this stationary 



block is therefore not affected by the approaching rupture 
front, and a spike will therefore appear at the rupture front. 
Spikes may also appear as an event arrests, also caused by 
one part of the interface slipping while another part is stuck. 

In all the following, the viscous damping introduced in 
this section will be used. 



3.2 Introducing a tangential stiffness of the interface 

As discussed above, the model by Maegawa et al. exhibits 
an unphysical scaling with N. It is possible to understand 
this scaling by considering how the system is tangentially 
loaded. As stated above, the driving force only acts on block 
1 . In order for an event to nucleate, this block has to reach its 
static friction threshold, which is proportional to the normal 
force: ji s pi °< p\ °= l/N. Since the added driving force per 
time is independent of N, the time between two events will 
be proportional to l/N, and the frequency of events con- 
sequently scales as N, an argument which is fully consis- 
tent with Fig. [5] The origin of the odd TV-dependence of the 
model by Maegawa et al. is therefore the decreasing size of 
the loading region as N is increased. 

In a physical system, like in the experiments by Mae- 
gawa et al, the loading region has a well defined spatial ex- 
tension, which is a combination of various effects. First, the 
tangential loading is applied at the trailing edge of the slider 
at some effective height h above the interface. As discussed 
experimentally in 1 1611 and modelled in [21], such a loading 
condition makes the tangential stress very high in the region 
near the trailing edge, the extension of which is of order h. 
Second, the interface between slider and base is not rigid. 
Both surfaces are rough and the multi-contact layer between 
them has a finite tangential stiffness associated to the tan- 
gential deformation of the microasperities involved in the 
contact. Such stiffness can be measured experimentally 83311 
and is found much smaller than the slider's bulk stiffness. 
Such a low interfacial stiffness can be responsible for devia- 
tions with respect to AC friction in static contacts [34]. The 
interfacial stiffness also results in a physical finite size of the 
loading region, since a localized tangential force at the inter- 
face will induce tangential strains within the rough layer not 
only at the loaded point but also in its neighbourhood. 

The first effect relates to the elastic coupling of points 
of the interface through the slider's bulk, which cannot be 
accounted for explicitly by a ID model like the one consid- 
ered here. In contrast, the effect of the interfacial stiffness 



can be introduced 113211 in a ID model in the following way. 
Each block is initially attached to the track by a spring with 
stiffness k t as seen in Fig. [9^. The spring connecting block 
n to the track has a breaking strength fl s p n (Fig. [9}}). When 
the spring to the track is detached the block is subject to the 
kinetic friction force ±/i^p„ (Fig. [9};). As the velocity of this 
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Fig. 9 Sketch of the behaviour of the spring between block n and the 
track, (a) The static friction force on block n is exerted through a spring 
of stiffness k t attached to the track, (b) As the block is moving, the 
spring is stretched until it reaches its breaking strength fl s p n . (c) When 
the spring is broken, block n is subject to the kinetic friction force flkPn ■ 
(d) As the block stops, the spring reattaches at x — such that the 
total force on the block is zero at the time of reattachment. 
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Fig. 10 Tangential force T normalised by the normal force p at initi- 
ation of the first precursor, when an interfacial stiffness is considered. 
The first block has just reached it threshold for slip (x/p = /!.,). The 
tangential force decays exponentially according to Eq. \A6\ . Results 
obtained using N = 100, 6 = and k, = 10 7 N/m. 



block reaches zero, the spring reattaches such that the total 
force on the block is at the time of reattachment (Fig.[9]l). 
The friction force is now given by 



fn = 



-k t [u n - u„) if attached, 
— sign (ii „) Hkp,, if detached, 



(9) 



where u^ ck is the attachment point of the spring between 
block « and the track. It is given by 



stick 



(10) 



where u„ and X n are the position of and tangential force on 
block n at the instant of its last reattachment to the track. 
This causes the total force on block n to be zero at the time 
of reattachment. The spring then detaches at the time step in 
which one finds 



-k t (u 



.stick 



> HsPn- 



(ID 



The equations of motion are still given by Eqs. (0. 

As the system is loaded tangentially, a finite region around 
the driving point is affected. This is illustrated in Fig. [10] 
where the tangential force at the time of nucleation of the 
first precursor is plotted. The length of this region depends 
on the stiffness k t of the springs between the blocks and the 
track. 

Assuming (i)N^$> 1, (ii) the length of the loading region 
to be much smaller than L and (iii) slow loading compared 
to the internal dynamics of the system, we calculate ana- 
lytically in App. [B] the tangential force profile at the time 
of nucleation of an event. No assumption is made for the 
tangential force profile at the time of arrest of the previous 
event, and the calculation is therefore valid for all events, 
not only for the first precursor shown in Fig. [10] The result 



is an exponential decay of the tangential force with a char- 
acteristic length Iq, which consequently is a measure of the 
size of the loading region. From App. [B] Iq is related to k t 
by 



ESL 



(12) 



In this expression, we recognize Nk t to be the total stiffness, 
k] ot , of the interface (N springs of individual stiffness k t in 
parallel), which is a measurable quantity in a given experi- 
mental setup. We then obtain the relation 



ESL/kf 



(13) 



which shows that the size Iq of the loading region is now 
independent of N, i.e. of the model's spatial resolution. 

In the model by Maegawa et al., the simple rigid-plastic- 
like AC friction ruled the behaviour of the interface. Now 
that the interfacial stiffness is introduced, the friction law 
is elasto-plastic-like, an improvement that has often been 
defended as a necessary extension of AC Ill9l.l33l - l3al . Our 
model is the simplest improvement of the model by Mae- 
gawa et al. that accounts for interfacial stiffness. 

The envelopes of the time dependence of both the tan- 
gential load and length of precursors in this improved model 
are not changed significantly compared to those obtained us- 
ing AC friction, as seen by comparing Figs. l2l and [TT1 How- 
ever, comparing Figs.fTT1tol4lshows a significant decrease in 
the number of events and a consequent increase in the am- 
plitude of the drops in Ft /f/v- The shear force profiles are 
also very similar to that of Fig. [8] except that spikes now 
decay exponentially with a characteristic length Iq. We will 
now check the model's behaviour with respect to its scaling 
with Af. To do this, one has to choose the interfacial stiffness 
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Fig. 11 Time evolution of Fj (a) and xr (b) in our improved model 
including both a relative viscous damping and an interfacial stiffness, 
using N = 100 and = 0. Comparison with Figs.|4]and|2]shows similar 
envelopes but very different numbers of events and amplitudes of the 
drops in /F^. 
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Fig. 12 Number of different kinds of events as a function of N, when 
the interfacial stiffness in taken into account. 8 = 0, f] = y/OAy/km and 
lo = 5 mm. (a) Number of global events for t 6 [5 s, 20 s] (solid line) 
and total number of events (dashed line), (b) Number of precursors 
longer than any previous one (solid line) and total number of precur- 
sors (dashed line). Comparison with Fig. \5\ shows that introduction of 
an interfacial tangential stiffness suppresses the resolution dependence 
of the numbers of global events and of precursors longer than any pre- 
vious one. 



kt (or equivalently a loading zone size Iq). The value of k t 
could be calibrated using an experimental measurement of 
the total interfacial stiffness £, tot . Here, we do not have ac- 
cess to such a measurement, so we rather exploit the fact 
that the number of precursors longer than any previous one 
is controlled by the parameter lo. By trial and error, we found 
that Iq = 5 mm, corresponding to k, = 10 7 N/m produced a 
number of such events similar to that observed in the exper- 
iments by Maegawa et al. . We will use this value of Iq in the 
following. 

In Fig. [T2]the number of events using the elasto-plastic 
friction law is plotted in the same way as in Fig. [5] for AC 
friction. Both the number of global events and the number 
of precursors longer than any previous one are now seen to 
be approximately constant. The total number of events and 
the total number of precursors are, however, still increasing 
with N. This is mainly due to events involving one block, i.e. 
events in which no real front propagation occurs. In other 
words, interfacial friction is found to satisfactorily solve the 
resolution dependence of the model's results, provided one 
considers events that have a measurable length. 



3.3 Introduction of an initial tangential force distribution 

The introduction of an internal viscosity and an interfacial 
stiffness in the model by Maegawa et al. allowed us to ob- 
tain force distributons and numbers of micro-slip events that 
were physically sound. However, the predicted L p vs. Ft 
curves still follow the same shapes as those shown in Fig. [7] 
which, as already mentioned, deviate significantly from those 
obtained experimentally. In an effort to further improve the 
ID model, we note that one of the main differences between 
the model and the experiment is the initial tangential force 



distribution (when no external tangential load has been yet 
applied). In the model, such forces are assumed to be zero 
all along the contact. However, in the experiments, both the 
slider and the base undergo different expansion rates dur- 
ing application of the normal loading. As discussed in 13711 . 
the associated slip at the interface is impeded by friction, 
thus yielding a significant tangential force distribution at the 
interface. Such distributions have been measured to be anti- 
symmetric 118], in agreement with basic contact mechanics 
calculations |6], thus ensuring Ft(0) = Ym=i T «(0) = 0. This 
effect is again a bulk effect which is quantitatively repro- 
duced in 2D models J2H1 . Here, in ID, we will only study 
the qualitative influence of an initial tangential force distri- 
bution on the length of precursors. We consider the simple 
linear distributions shown in Fig. IT3h . 

Implementation of an initial tangential force requires an 
initial relative displacement of the blocks. The initial tan- 
gential forces are given by 

f Jt(M2(0)-Mi(0))+Fr(0), «=1 
T„(0)= I k(u n+1 (0)-2u n (0) + u n ^(0)),2<n<N-l 
[*(it iv , -i(0)-Mjv(0)), n = N, 

(14) 

and by choosing Ft(0) = and u\ (0) = the above equation 
can be rewritten to 

fO, n = l 

k„(0) = I ui+Xn/k, n = 2 (15) 

[ 2u„-i-u„-2 + %-i/k, n = 3,4,..., N, 

thus enabling calculation of the initial positions of all blocks 
given T„. The initial attachment position of the interfacial 
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Fig. 13 (a) Three different antisymmetric initial tangential force distri- 
butions aimed at modelling the effect of friction-frustrated differential 
Poisson expansion during normal loading, o: homogeneous distribu- 
tion, x and +: linear profiles of increasing slope (/3 = 0.225 and 0.45 
in Eq. ( 1521 respectively), (b) The length of precursors corresponding to 
the three initial tangential force profiles shown in (a). Results obtained 
using N = 100, 6 = 0, T] = \/0.1\/km and Iq — 5 mm. Solid lines are 
the analytical predictions of the precursor lengths discussed in Sec. [4] 



springs is then calculated using Eq. (flOl) . which ensures that 
the total force on each block is zero at t = 0. 

The length of precursors corresponding to the three ini- 
tial tangential force profiles shown in Fig. [T3h is shown in 
Fig.[T3b. using = 0. With respect to a zero initial force dis- 
tribution, the stronger the asymmetry the lower the L p vs. Fj 
curve, i.e. the shorter the precursors for the same tangential 
force. Moreover, the slope of the curves is significantly in- 
creased at large forces. Both effects lead to a significantly 
improved qualitative agreement with the experimental re- 
sults by both Maegawa et al. and Rubinstein et al. QjJ.This 
shows that the initial force distribution, which arises natu- 
rally in 2D models 112 ill , is a crucial parameter for the kine- 
matics of precursors to sliding. Note that the number of pre- 
cursors, which is closely related to the choice of Iq (as men- 
tioned in the previous subsection) is only weakly affected by 
the introduction of an initial tangential force distribution. 



4 Analytical prediction of precursor lengths 

In order to complete this study we derive an analytic predic- 
tion for the precursor length as a function of the tangential 
force at event arrest, L p vs. Fj. We will first look at the sim- 
pler case of the model by Maegawa et al., and then extend 
the prediction to our improved model using both an interfa- 
cial stiffness and an initial tangential force profile. 



4. 1 Prediction in the model by Maegawa et al. 

Assume that a precursor has reached block n p and has the 
length L p = (n p /N)L. We want to calculate the tangential 
force Ft at the time of arrest of this event. At that time, all 
blocks are stuck, so that 



Fj 



N 
n=\ 



(16) 



This means that, given the tangential force distribution at 
event arrest, the corresponding tangential force is found us- 
ing Eq. ( fT6b . 

According to Figs. [6] and [8] the tangential force is ob- 
served to be approximately equal to the kinetic friction level 
from block 1 to n p , and elsewhere. Using this assumption, 
Eq. ( TTol i yields 



f t = £ Pn 



2n-N-l 
N- 1 



-e. 



(17) 



/;=1 n=l 

where Eq. (f5J) has been inserted for p n . If N 3> 1 , the sum can 
be approximated by an integral and n replaced by x = nL/N, 
which yields 



Fj ■ 



F N N 
'^NLJo 



1 



2(xN/L)-N- 1 



N-l 



dx, 



and approximating N±l ~N yields 

Fn r L p 

F T KH k -p- 1- 2 */L -1 dx 
L Jo 



Ft ~ Mit-fiV 



L 



1 



1 -■ 



(18) 

(19) 
(20) 



As seen in Fig. [7] this prediction is in very good agree- 
ment with our simulation results. The deviation between the 
actual precursors and the analytical curve is the result of a 
slightly incorrect assumed tangential force profile. Note that 
a similar good prediction scheme, numerical rather than ana- 
lytical, was previously developed in the 2D study by Tr0mborg 
et al. lEl]. 

According to Eq. ( f20b the global static friction coeffi- 



cient n s (jj, s 



is easily evaluated for L p = L from Eq. (1201 ) 



is independent of 9 and almost equal to the local kinetic 
friction coefficient This is in agreement with the model 
result in Fig. [7] and with the results of the 2D model by 
Tr0mborg et al. l2lll . 
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4.2 Prediction in our improved model 

Prediction of the precursor length in our improved model 
including both an interfacial stiffness and an initial tangen- 
tial force profile follows the same line as that for the model 
by Maegawa et al. However, blocks may now move even 
though no event is occurring. Despite this, it is expected that 
block accelerations are small when all track springs are at- 
tached, which leads to the approximate validity of Eq. ( TToT ). 

An approximate tangential force profile at the time of 
arrest of an event (Fig. fl4l has to be found for a given L p . 
Again, blocks in [0,L p ] are assumed to have a tangential 
force equal to the kinetic friction level. Blocks in the in- 
terval [L P ,L], however, now needs to be taken into account 
for two reasons: both the initial tangential force profile and 
the springs to the track lead to a non-zero tangential force 
for x > L p . The form of this profile in a static situation has 
been calculated in App. [B] and is given in Eq. d46l i. where 
T°(x) now is the initial tangential force profile. However, it 
has to be modified to take into account that the loaded block 
is not located at x = 0, but at x = L p , and that this block does 
not have to be loaded up to the static friction level, but may 
take some other value, say (Xp„ , with a a coefficient to be 
defined. Our assumed tangential force profile at the arrest of 
a precursor of length L p is therefore given by 



( Hkp(x), 
r(x) = < 

I (ccp(L p )- 



x e [o,l p ] 

T°(x))e~^T +T°(x), x g [L pi L] 



(21) 



We have considered that all values of the amplitude of the 
peak at x = L p have the same probability to occur between 
pLkP and jj, s p, so that we have given a its average value 
a = (jj. s + jUjfc)/2 in the predictions seen in Fig. [T3b . Fig- 
ure [14] shows the assumed tangential force profile and the 
actual tangential force profile at the arrest of an event, and 
the agreement is seen to be satisfactory. 

From Eqs. (f2TT >. one can calculate Fj{L p ) in the same 
way as was done for the model by Maegawa et al. This cal- 
culation is provided in AppendixICl As seen in Fig.fTjb. the 
prediction scheme works well. Deviations between the ac- 
tual precursors and the analytical curve have two contribu- 
tions: incorrectly assumed tangential force profile and iner- 
tial effects, where the former gives the largest contribution. 
Again, the global static friction coefficient jj.s is seen to be 
approximately equal to the local kinetic friction coefficient 
jU£. The prediction curves appear to bend slightly backwards 
at Fj/Fn ~ 0.45. In our continuous prediction scheme, this 
corresponds to micro-slip fronts that are so close to the lead- 
ing edge that all blocks at x > L p have a tangential force 
above the kinetic friction level, causing Ft to get smaller 
as the front moves further. In reality, however, no precursor 
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Fig. 14 The tangential force profile at the arrest of the 14th precursor 
plotted in Fig.[T3](red dots) and the predicted profile used in the analyt- 
ical prediction of precursor lengths (black drawn line) using N = 100, 
9 = 0,1] = VOAVkm, / = 5 mm and j8 = 0.225 in Eq. {55]l. 

ever arrests in such a state but propagates all the way to the 
trailing edge. 



5 Conclusion 

Recent experimental results about the transition from static 
to kinetic friction in line contacts have triggered the study 
of the deterministic dynamics of ID spring -block friction 
models in which driving is applied at one extremity of the 
chain of blocks. In this Letter, we have improved the sim- 
plest of such models 111711 in order to solve its intrinsic un- 
physical resolution-dependence and to ameliorate its quali- 
tative agreement with experimental results on the kinematics 
of micro-slip fronts along the contact. In particular, the intro- 
duction of a tangential stiffness of the interface, by introduc- 
ing a new length scale in the model, practically suppresses 
its resolution dependence and allows for reproduction of re- 
alistic numbers of precursory micro-slip fronts. The addi- 
tional introduction of an initial tangential force distribution 
at the interface significantly improves the agreement with 
the evolution of the precursor length with the external tan- 
gential load obtained in experiments. Our improved model 
is intented to serve as a framework for more complex fric- 
tion models when robust comparisons with experiments are 
desired. 

We focused on ID models because they are simple enough 
to enable deep insights into the qualitative effects of the 
model's parameters. However, it is known that 2D models 112 111 
are required to provide quantitative agreement with experi- 
ments. In this respect, the improvements brought to the ID 
model are effective ways to account for intrinsically 2D ef- 
fects: First, the length scale introduced through the interfa- 
cial stiffness enables coupling between remote points along 
the interface, analogous to the coupling through the slider's 
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bulk; Second, the initial tangential force distribution accounts 
for the shear stress arising form the differential Poisson ex- 
pansion of two bodies pressed together. 

In analogy with 2D results we developed, based on the 
well defined force distribution left by an arrested precursor, 
an efficient analytical prediction for the precursors' length 
as a function of the external tangential load applied. We also 
find that, like in 2D, the macroscopic static friction coeffi- 
cient of a side-driven contact is approximately equal to its 
microscopic kinematic friction coefficient. 
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A Relative viscous damping in a linear chain of blocks 

If friction forces are ignored, the equation of motion for an infinite 
chain of blocks connected by springs is given by 



mii„ =k(u n+ i -2m„ + m„_i) + J7(»„ + i -2k„ + 
We then assume a solution of the form 

u n (t)=e^'e iKna , 



(22) 



(23) 



where f K e C and f gl. Inserting Eq. i23\ into Eq. i22\ yields the 
relation 



m£ = k (e iKa -2 + e-' Ka ) + 77 £ K (e iKa -2 + e - iK ") 



which can be simplified to 

m£ + 47] sin 2 — J £ K + 4k sin 2 ( — 



(24) 



(25) 



B Tangential force profiles and characteristic length 
with a tangential stiffness of the interface 

An analytical expression for the characteristic length Iq can be found. 
In order to do so, the following assumptions are made: N^> 1,/o/LCl 
and slow loading compared to the internal dynamics of the system, 
which enables a static analysis. The system is first placed in a static 
state with an initial shear force profile given by T„, and then loaded 
slowly from the left. The equilibrium of all non-edge blocks writes 



* (w n+ i - 2u„ +u n - l )-k, (u n - uf ck ) = 0. 
We introduce a new variable u' n defined by 



u„ = u n + u n , 



(30) 



(3D 



where u„ is the initial position of block n. Inserting Eq. d3U into Eq. d30t 
yields 



k « +1 - 2u' n + u' n _ , ) - k,u' n + f° - k, (u° n - uf A ) = 0, 



where 

T° = *(«° +1 -2 M ° + ul_ x ). 



(32) 



(33) 



The two terms T° and —k, (wj] — «, s 1 " ck ) cancel in Eq. <32b since the 
initial state is static, and thus 



k(u' n+l -2u'„+u' n _ x )-k,u' n =Q. 
The above equation can be rewritten to 

ka 2 < +l -H + <- l _ ktU > n = 0i 



(34) 



(35) 



where a = L/(N X — 1) is the lattice spacing. Since N 1, the first 
term in Eq. i35\ can be replaced by the second spatial derivative, and 
replacing u'„ with u' (na) =u! (x) yields 



(36) 



' = -4 sin 2 



The complex parameter £ K is then given by 



-4rj sin 2 (f ) ± y / 16rj 2 sin 4 (f^) - 16/tmsin 2 (§^) 
2m 



(27) 



The system is critically damped when Eq. ((25} only has one solution 
for £,e, which occurs when the square root is zero: 



T) sin 



: km 



kin 



(28) 



The oscillations that are to be reduced have a wavelength X = 2a, i.e. 
a wave number K = 2lt/X = it/a. Inserting this into Eq. )28b leads to 



r\ c = Vkm, 



(29) 



which is the value of the damping coefficient 77 for which waves of 
wavelength X = 2a are critically damped. Since the absolute value of 
sin is always smaller than one, choosing 77 = \f~km will cause all other 
waves to be under-damped. 



(26) which has the general solution 



u'(x) ^Ae*!' +Be- xl '°, l = \l -a. 

V k, 

The shear force is given by 

= k («,',+ 1 " 2",', + ) + t°. 
By replacing again finite differences with second order derivatives, 

, •> d 2 u'(x) n , , 
z(x) =ka 2 I ' +t°(x), 



(37) 



(38) 
(39) 



(40) 



and the general expression for the shear force profile can be found by 
using Eq. (37) , which yields 



(41) 



The system is loaded from the left, and at the beginning of an event 
the shear force on block 1 is equal to the static friction threshold fi s p\ . 
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Provided Iq/L<S. 1, the trailing edge will not be affected by the loading. 
The latter of these two boundary conditions yields 



t(L) = ^ (a^ +Be- L "«) + T°(L) 



ka l ; 

Ad 



2 i 2 / \ 



i.e. A = 0. The first boundary condition yields 



(42) 

(43) 
(44) 



The above integrals can be calculated easily, and the result is the tan- 
gential load Fj as a function of the precursor length L p : 



F T (L„)=F N 
/o 



L 

p(l- 



L- 

L—L, 



L 



We observe that again Fj (L) = fi^F^ 



1 +, 



. {L-Lp)L p 
V- 



L-Lj, 

- a I 1 — e *o 



(54) 



t(0) = ^5 + t°(L) = ^i, 

and the shear force is therefore given by 

z{x) = (n sPi - T°(x)) e-'I'o + T°(x). 



The characteristic length lg is given by Eq. 437b . and inserting for 
k given by Eq. j4j and a yields 



k 

—a = 



ELS 



k t \j(N-l)kr 
and hence Eq. )12t for N > > 1 . 



(47) 



Note that in a 3D situation, the exponential decay of the tangential 



stress with x would be replaced by a power law 



decay or t 

[US. 



C Derivation of the prediction of precursor lengths in 
our improved model 

We start with Eq. U6b and use the assumed shear force profile in Eq. 12 It . 
shown in Fig.[l4] Again we go to the limit N — > °°, resulting in the sub- 
stitution 



N 

E t « 

n=l 



!\{x) 

Jo 



N i-l 

T 



dx. 



■ xN/L. 



(48) 



The tangential force after a precursor of length L p is then given by 



N 

Fr = _ 
_ N 



t(x)dx+ I t(x)cLc 



L " H k p(x)dx+ [ L (ap(L p ) - t°(x)) e •<>" +T°(x)dx 
Jo Jl„ 



(50) 

We limit ourselves to predicting the precursors in Fig. 1 13b . i.e. using 
a tangential interfacial stiffness and a linear initial shear forces as de- 
picted in Fig. 113b . but with 0=0. The normal and initial shear force 
are then given by 



p(x) = p = Fjy / N = constant, 



2{x-L/2) 



where the parameter /3 determines the slope in the initial shear force 
profile. Inserting Eqs. <5U and J521 into Eq. l !50t yields 



Fj : 



N 



H k pdx+ 



LX ap - pp 



2(x-L/2)\ -'J±t 2(x-L/2) j 
1 ' x e ^T+jlp— ^-dx 
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